///@file Polynomial.cpp
///@brief Implementation of the class representing a Polynomial
///@author Arnaud Duval
///@version 0.0
///@date 2009/02/04

#include "Polynomial.h"
#include <cstring>
#include <cmath>
#include "../../utils/include/StringUtils.h"


///@brief Default constructor
Polynomial::Polynomial()
{
	degree = -1;
}

///@brief Constructor precising the degree of the polynomial
///@param d Degree of the polynomial
Polynomial::Polynomial(int d)
{
	assert(d >= 0);
	degree = d;
	a = new double[d+1];
	for(int i = 0 ; i <= d ; i++)
		a[i] = 0.;
}

///@brief Constructor of specific polynomial
///@param d Degree of polynomial
///@param t Type of polynomial : 'H' for Hermite, 'L' for Legendre
Polynomial::Polynomial(int d, char t)
{
	assert(d >= 0);
	degree = d;
	a = new double[d+1];
	if(t == 'H')	//Construction of Hermite Polynomial
	{
		Polynomial *serie;
		serie = new Polynomial[d+1];
		serie[0] = Polynomial(0);
		(serie[0])[0] = 1.;
		if(degree > 0)
		{
			serie[1] = Polynomial(1);
			(serie[1])[0] = 0.;
			(serie[1])[1] = 1.;
		}
		if(degree > 1)
		{
			for(unsigned int i = 2 ; i <= degree ; i++)
				serie[i] = serie[1]*serie[i-1] - ((double)(i-1))*serie[i-2];
		}

		for(unsigned int i = 0 ; i <= degree ; i++)
			a[i] = (serie[degree]).a[i];
		delete[] serie;
	}
	else if(t == 'L')
	{
		Polynomial *serie;
		serie = new Polynomial[d+1];
		serie[0] = Polynomial(0);
		(serie[0])[0] = 1.;
		if(degree > 0)
		{
			serie[1] = Polynomial(1);
			(serie[1])[0] = 0.;
			(serie[1])[1] = 1.;
		}
		if(degree > 1)
		{
			for(unsigned int i = 2 ; i <= degree ; i++)
				serie[i] = ((2.*(i-1)+1.)*serie[1]*serie[i-1] - (((double)(i-1))*serie[i-2]))*(1./((double)(i)));
		}

		for(unsigned int i = 0 ; i <= degree ; i++)
			a[i] = (serie[degree]).a[i];
		delete[] serie;
	}
}

///@brief Copy constructor
///@param p Polynomial to copy
Polynomial::Polynomial(const Polynomial &p)
{
	degree = p.degree;
	a = new double[degree+1];
	for(unsigned int i = 0 ; i <= degree ; i++)
		a[i] = p.a[i];
}

///@brief Copy constructor for temporary object
///@param p Polynomial to copy
Polynomial::Polynomial(Polynomial* p)
{
	degree = (*p).degree;
	a = new double[degree+1];
	for(unsigned int i = 0 ; i <= degree ; i++)
		a[i] = (*p).a[i];
}

///@brief Destructor
Polynomial::~Polynomial()
{
	if(degree >= 0)
		delete[] a;
}

///@brief Overload of = operator
///@param p right hand operand
Polynomial& Polynomial::operator=(const Polynomial &p)
{
	if(degree >= 0)
		delete[] a;
	degree = p.degree;
	a = new double[degree+1];
	for(unsigned int i = 0 ; i <= degree ; i++)
		a[i] = p.a[i];
	return (*this);
}

///@brief Overload of = operator for temporary object
///@param p right hand operand
Polynomial& Polynomial::operator=(Polynomial* p)
{
	degree = (*p).degree;
	a = new double[degree + 1];
	for(unsigned int i = 0 ; i <= degree ; i++)
		a[i] = (*p).a[i];
	return (*this);
}
///@brief Overload of &= operator : copy without size check
///@param p right hand operand
Polynomial& Polynomial::operator&=(const Polynomial &p)
{
	if(p.degree < 0)
	{
		if(degree >= 0)
			delete[] a;
		degree = -1;
	}
	else
	{
		if(p.degree != degree)
		{
			if(degree >= 0)
				delete[] a;
			degree = p.degree;
			a = new double[degree+1];
		}
		memcpy(a, p.a, (p.degree+1)*sizeof(double));
	}
	return (*this);
}

///@brief Overload of [] operator
///@param i index of returned coefficient
double& Polynomial::operator[](unsigned int i) const
{
	assert((i >= 0)&&(i <= degree));
	return a[i];
}

///@brief Evaluate polynomial value
///@param x Value of x for while the polynomial must be evaluated
double Polynomial::Evaluate(double x)
{
	double res = 0;
	for(unsigned int i = 0 ; i <= degree ; i++)
		res += a[i] * pow(x, (double)i);
	return res;
}

///@brief Return string with polynomial expression
///@param var Expression of polynomial variable
std::string Polynomial::Expression(std::string var)
{
	std::string res;
	assert(degree >= 0);
	res = DoubleToString(a[0], 5);
	for(unsigned int i = 1 ; i <= degree ; i++)
		res += " + "+DoubleToString(a[i], 5)+" "+var+"^"+IntToString(i, 3);
	return res;
}

///@brief Overload of * operator : product polynomial-scalar
///@param p left hand operand : polynomial
///@param s right hand operand : scalar
Polynomial operator * (const Polynomial& p, const double& s)
{
	assert(p.degree >= 0);
	Polynomial res(p.degree);
	for(unsigned int i = 0 ; i <= p.degree ; i++)
		res.a[i] = p.a[i] * s;
	return res;
}

///@brief Overload of * operator : product scalar-polynomial
///@param s left hand operand : scalar
///@param p right hand operand : polynomial
Polynomial operator * (const double& s ,const Polynomial& p)
{
	assert(p.degree >= 0);
	Polynomial res(p.degree);
	for(unsigned int i = 0 ; i <= p.degree ; i++)
		res.a[i] = p.a[i] * s;
	return res;
}
///@brief Overload of * operator : product polynomial-polynomial
///@param p1 left hand operand : polynomial
///@param p2 right hand operand : polynomial
Polynomial operator * (const Polynomial& p1, const Polynomial& p2)
{
	assert(p1.degree >= 0);
	assert(p2.degree >= 0);
	Polynomial res(p1.degree + p2.degree);
	for(unsigned int i = 0 ; i <= p1.degree ; i++)
	{
		for(unsigned int j = 0 ; j <= p2.degree ; j++)
			res.a[i+j] += p1.a[i] * p2.a[j];
	}
	return res;
}

///@brief Overload of + operator
///@param p1 left hand operand
///@param p2 right hand operand
Polynomial operator + (const Polynomial& p1, const Polynomial& p2)
{
	assert(p1.degree >= 0);
	assert(p2.degree >= 0);
	Polynomial res;
	if(p1.degree > p2.degree)
		res = Polynomial(p1.degree);
	else
		res = Polynomial(p2.degree);
	for(unsigned int i = 0 ; i <= res.degree ; i++)
	{
		if(((p1.degree >= p2.degree)&&(i <= p2.degree))||((p2.degree >= p1.degree)&&(i <= p1.degree)))
			res.a[i] = p1.a[i] + p2.a[i];
		else
		{
			if(p1.degree > p2.degree)
			{
				res.a[i] = p1.a[i];
			}
			else
			{
				res.a[i] = p2.a[i];
			}
		} 
	}
	return res;
}


///@brief Overload of - operator
///@param p1 left hand operand
///@param p2 right hand operand 
Polynomial operator - (const Polynomial& p1, const Polynomial& p2)
{
	assert(p1.degree >= 0);
	assert(p2.degree >= 0);
	Polynomial res;
	if(p1.degree > p2.degree)
		res = Polynomial(p1.degree);
	else
		res = Polynomial(p2.degree);
	for(unsigned int i = 0 ; i <= res.degree ; i++)
	{
		if(((p1.degree >= p2.degree)&&(i <= p2.degree))||((p2.degree >= p1.degree)&&(i <= p1.degree)))
			res.a[i] = p1.a[i] - p2.a[i];
		else
		{
			if(p1.degree > p2.degree)
			{
				res.a[i] = p1.a[i];
			}
			else
			{
				res.a[i] = -p2.a[i];
			}
		} 
	}
	return res;
}

///@brief Overload of << operator
///@param s out stream
///@param p Polynomial object to put in the out stream
std::ostream& operator<<(std::ostream &s, const Polynomial& p)
{
	assert(p.degree >= 0);
	s << p.a[0];
	for(unsigned int i = 1 ; i <= p.degree ; i++)
		s << " + " << p.a[i] << " X^" << i;
	return s;
}